#!/usr/bin/env python
"""

    single_linkage_cluster.py
    [--log_file PATH]
    [--verbose]

"""

################################################################################
#
#   single_linkage_cluster
#
#
#   Copyright (c) 7/5/2010 Leo Goodstadt
#
#   Permission is hereby granted, free of charge, to any person obtaining a copy
#   of this software and associated documentation files (the "Software"), to deal
#   in the Software without restriction, including without limitation the rights
#   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#   copies of the Software, and to permit persons to whom the Software is
#   furnished to do so, subject to the following conditions:
#
#   The above copyright notice and this permission notice shall be included in
#   all copies or substantial portions of the Software.
#
#   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
#   THE SOFTWARE.
#################################################################################

import sys, os
from itertools import izip


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Functions


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#_________________________________________________________________________________________

#   to_histogram

#_________________________________________________________________________________________

def nth_percentiles (array, N = 5):

    if len(array) == 0:
        array[:] = [0] * (N + 1)
        return
    s_array = sorted(array)
    sz = len(array)
    reduced_values = []
    for i in range(21):
        pos = min(sz - 1, int(sz * i * N * 0.01))
        reduced_values.append(s_array[pos])
    array[:] = reduced_values


import logging
from quantile_histogram import quantile_histogram
def analyse_cluster (clusters, logger, logger_level = logging.INFO):
    lengths = sorted(len(c) for c in clusters)
    logger.log(logger_level, "%8d clusters" % len(clusters)                          )
    logger.log(logger_level, "%8d items   " % sum(lengths)                           )
    logger.log(logger_level, "%8d items in largest cluster" % max(lengths)         )
    logger.log(logger_level, "%8d items median" % lengths[len(clusters) / 2] )
    logger.log(logger_level, "cluster size histogram")
    quartile_sz = 100.0 / 10
    for i, h in enumerate(quantile_histogram (lengths, 10)):
        logger.log(logger_level, "  %5d%% -- %10d" % (i * quartile_sz, h))


#
#   expects list of identifiers
#   list of identifier pairs which are linked (only first two fields are retained)
#
#
def __lookup_name (name, names, names_to_index, clusters, id_to_cluster_id):
    """
    If name in names
        return the id and cluster_id for name
    Else:
        Add to list of names and as a member of new cluster
        and return the id and cluster_id for name
    """

    if not name in  names_to_index:
        #
        #   add to list of names
        #
        name_id = len(names)
        names.append(name)
        names_to_index[name] = name_id

        #
        #   add as only member to new cluster
        #
        cluster_id = len(clusters)
        clusters.append([name_id])
        id_to_cluster_id[name_id] = cluster_id
    else:
        name_id = names_to_index[name]
        cluster_id = id_to_cluster_id[name_id]

    return name_id, cluster_id

def single_linkage_cluster (name_tuples):

    #
    #   list of names
    #
    names    = []
    names_to_index = dict()

    #
    #   start with each orthoset in its own cluster
    #
    clusters            = []
    id_to_cluster_id    = dict()



    for name_tuple in name_tuples:

        #   convert links to indices
        name = name_tuple[0]

        #
        # This is the cluster where everything will end up
        #
        id0, cluster_id0 = __lookup_name (name, names, names_to_index, clusters, id_to_cluster_id)

        #print "destination %d in cluster %d" % (id0, cluster_id0)

        for other_name in name_tuple[1:]:
            other_id, other_cluster_id = __lookup_name (other_name, names, names_to_index, clusters, id_to_cluster_id)

            # part of same cluster already: ignore
            if other_cluster_id == cluster_id0:
                continue

            #
            # merge contents of other_cluster_id with cluster_id0
            #
            if clusters[other_cluster_id] == None:
                raise Exception ("cluster #%d for id %d is none" % (other_cluster_id, other_id))
            if clusters[cluster_id0] == None:
                raise Exception ("Destination cluster #%d for id %d is none" % (cluster_id0, id0))
            for other_cluster_member in clusters[other_cluster_id]:
                id_to_cluster_id[other_cluster_member] = cluster_id0
                #print "  Move id %d in cluster %d -> %d " % (other_cluster_member, other_cluster_id, cluster_id0)

            # move contents of other_cluster_id to cluster_id0;
            clusters[cluster_id0].extend(clusters[other_cluster_id])
            #print "  cluster %d now has %s " % (other_cluster_id, clusters[other_cluster_id])
            clusters[other_cluster_id] = None

    # translate ids back to names
    #print clusters
    return sorted(sorted(names[id] for id in c) for c in clusters if c != None)


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Main logic


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
if __name__ == '__main__':
    import logging,sys
    import logging.handlers

    MESSAGE = 15
    logging.addLevelName(MESSAGE, "MESSAGE")

    def setup_std_logging (logger, log_file, verbose):
        """
        set up logging using programme options
        """
        class debug_filter(logging.Filter):
            """
            Ignore INFO mesages
            """
            def filter(self, record):
                return logging.INFO != record.levelno

        class NullHandler(logging.Handler):
            """
            for when there is no logging
            """
            def emit(self, record):
                pass

        # We are interesting in all messages
        logger.setLevel(logging.DEBUG)
        has_handler = False

        # log to file if that is specified
        if log_file:
            handler = logging.FileHandler(log_file, delay=False)
            handler.setFormatter(logging.Formatter("%(asctime)s - %(name)s - %(levelname)6s - %(message)s"))
            handler.setLevel(MESSAGE)
            logger.addHandler(handler)
            has_handler = True

        # log to stderr if verbose
        if verbose:
            stderrhandler = logging.StreamHandler(sys.stderr)
            stderrhandler.setFormatter(logging.Formatter("    %(message)s"))
            stderrhandler.setLevel(logging.DEBUG)
            if log_file:
                stderrhandler.addFilter(debug_filter())
            logger.addHandler(stderrhandler)
            has_handler = True

        # no logging
        if not has_handler:
            logger.addHandler(NullHandler())


    #
    #   set up log
    #
    logger = logging.getLogger("me")
    setup_std_logging(logger, "me.log", True)
#   debug code not run if called as a module
#
    #88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

    #   Testing


    #88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

    import unittest
    from random import randint, shuffle
    class Test_single_linkage_cluster(unittest.TestCase):

        #       self.assertEqual(self.seq, range(10))
        #       self.assert_(element in self.seq)
        #       self.assertRaises(ValueError, random.sample, self.seq, 20)
        #
        @staticmethod
        def create_clusters (cnt_names, cnt_clusters):
            import random
            clusters      = [[] for i in range(cnt_clusters)]

            names = ["name_%03d" % i for i in range(cnt_names)]

            # sort into clusters
            for n in names:
                dest_cluster = cnt_clusters
                while dest_cluster >= cnt_clusters:
                    dest_cluster = int(cnt_clusters * random.lognormvariate(1, 5) * 0.5)
                    #dest_cluster = int(random.lognormvariate(cnt_clusters / 2, 5))
                clusters[dest_cluster].append(n)

            # ignore empty clusters and sort
            clusters = sorted(sorted(c) for c in clusters if len(c) > 1)
            analyse_cluster (clusters, logger, MESSAGE)
            return clusters, names

        def test_pairwise(self):
            """
                test pairwise tuples
            """
            print >>sys.stderr, "Setup..."
            cnt_sim_names    = 20000
            cnt_sim_clusters = 3000
            cnt_sim_tuples   = 100000
            #cnt_sim_names    = 20
            #cnt_sim_clusters = 4
            #cnt_sim_tuples   = 8
            clusters, names = self.create_clusters(cnt_sim_names, cnt_sim_clusters)

            # assign links
            name_tuples = []


            # add random links
            cnt_extra_tuples = cnt_sim_tuples
            while len(name_tuples) < cnt_extra_tuples:
                cluster_id = randint(0, len(clusters) - 1)
                cluster = clusters[cluster_id]
                if len(cluster) == 1:
                    continue

                target_tuple_size = randint(2, len(cluster))
                #target_tuple_size = 2

                name_tuple = set()
                while len(name_tuple) < target_tuple_size:
                    name_id  = randint(0, len(cluster) - 1)
                    name_tuple.add(cluster[name_id])

                name_tuples.append(list(name_tuple))

            # add pairwise to ensure coverage
            for cluster in clusters:
                if len(cluster) == 1:
                    continue

                # pairwise between successive cluster members too ensure that we can recover cluster
                for i in range(len(cluster) - 1):
                    name_tuples.append((cluster[i], cluster[i+1]))

            shuffle(name_tuples)


            print >> sys.stderr, "setup done: SLC..."
            import time
            start_time = time.time()
            relinked_cluster = single_linkage_cluster (name_tuples)
            print >>sys.stderr, "SLC done in %.2fs" % (time.time() - start_time)
            if relinked_cluster != clusters:
                print  >>sys.stderr, relinked_cluster
                print  >>sys.stderr, clusters
            #for t in name_tuples:
            #    print t
            self.assertEqual(relinked_cluster, clusters)


    #
    #   call unit test without parameters
    #

    if sys.argv.count("--debug"):
        sys.argv.remove("--debug")
    unittest.main()









